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We have used quantum Monte Carlo methods to calculate the zero-temperature phase diagram of 
the two-dimensional homogeneous electron gas. We find a transition from a paramagnetic fluid to an 
antiferromagnetic triangular Wigner crystal at density parameter = 31(1) a.u. and a transition to 
a ferromagnetic crystal at = 38(5) a.u. The fully spin-polarized fluid is never stable. We searched 
for, but did not flnd, the ferromagnetic "hybrid" phase proposed by H. Falakshahi and X. Waintal 
[Phys. Rev. Lett. 94, 046801 (2005)]. 
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The two-dimensional (2D) homogeneous electron gas (HEG) is one of the most important model systems in 
condensed-matter physics. It is our basic model of the conduction electrons in layered semiconductor devices, and 
consists of a set of electrons moving in 2D in a uniform, inert, neutralizing background. At high densities the HEG 
exists in a Fermi fluid phase, but at low densities it forms a Wigner crystal to minimize the electrostatic repulsion 
between the electrons l| . Classical 2D Wigner crystals can be produced in the laboratory by spraying electrons onto 



droplets of liquid He 



2l , and quantum Wigner crystals can be formed at the interface between two semiconductors 



As well as being objects of fundamental scientific interest, 2D Wigner crystals may be of use in quantum computing 
devices |4|. 

At high densities the fluid ground state is paramagnetic, but at low densities this phase is unstable to a ferromagnetic 
fluid. There is some experimental evidence that a ferromagnetic fluid is stable at densities between the paramagnetic 
fluid and crystal phases [5[, but these data are not conclusive. Previous theoretical studies [6|, |7| have reported a 
region of stability for the ferromagnetic fluid. However, the energy differences between the various phases are small, 
and highly accurate calculations are required to resolve them. Another issue of recent interest is the possible existence 



SHliJ. In this paper we address the issues of 



of intermediate phases in the vicinity of the fluid-crystal transition 
the stability of the ferromagnetic fluid, the fluid-crystal transition, the magnetic behavior of the crystal, and the 
possibility of a "hybrid" phase of the type proposed in Ref. 11 1. 



We have performed quantum Monte Carlo (QMC) 



12 



13| calculations for the 2D fluid and Wigner crystal phases, 



achieving much smaller statistical error bars than previous studies y, [7 



141. We used the variational and diffusion 



quantum Monte Carlo (VMC and DMC) methods as implemented in the CASINO code 



15j . DMC is the most accurate 



method available for studying quantum many-body systems such as electron gases. Fermionic symmetry is imposed 



via the fixed- node approximation [16[ , in which the nodal surface is constrained to equal that of a trial wave function 



The Ewald method was used to calculate the electron-electron interaction energy 



17[. We used extremely flexible 



Slater- Jastrow-backflow [18| wave functions, and wc optimized the crystal orbitals by directly minimizing the DMC 
energy. Our fixed-node DMC energies are therefore more accurate than those of earlier calculations. Finally, wc have 



19|. 



dealt with finite-size effects in a more sophisticated fashion than previous studies 

It is well-established that the triangular crystal has a considerably lower energy than competing lattices in the 
density range of interest. A triangular lattice was therefore used in all our crystal calculations. Antiferromagnetism 



is frustrated on a triangular lattice, and antiferromagnctic crystals are expected to form a spin liquid 20|. Instead, 
we have considered antiferromagnctic crystals consisting of alternating lines of spin-up and spin-down electrons. 



We used the Jastrow factor proposed in Ref. [2l|. Plane- wave orbitals exp(iG • r) were used in the Slater wave 
function in our fiuid calculations, while identical Gaussian orbitals exp(— Cr^) centered on lattice sites were used in 
our crystal calculations. The Slater wave function was evaluated at quasiparticle coordinates defined by the backflow 
(BF) transformation described in Ref. 



within VMC by minimizing the energy 



18|. The free parameters in the Jastrow factor and BF function were optimized 
22|. The BF functions for antiparallcl spins arc much larger than those for 
parallel spins, which arc already kept apart by Pauli exclusion. BF is less important in the crystal because the 
localization of the orbitals keeps the electrons apart. For example, at Vs ~ 30 a.u. [23|, BF lowers the DMC energy 



of the 42-electron paramagnetic fluid by 0.000036(3) a.u. per electron, whereas it lowers the energy of the 64-electron 
ferromagnetic Wigner crystal by only 0.0000023(3) a.u. per electron. Time-step bias was removed from our final DMC 
energies by linear extrapolation to zero time step 



24 



A target population of at least 1500 configurations was used 



24| 



in our production runs, making population-control bias negligible 

We optimized the Gaussian exponent C in the crystal orbitals by minimizing the Slater- Jastrow DMC energy. 
At Ts = 40 a.u., the DMC-optimized exponent of the ferromagnetic crystal (C^j^q w 0.0003 a.u.) is smaller than 
the exponent obtained by minimizing the VMC energy {Cyf^Q ~ 0.0006 a.u. with our Jastrow factor), which is, 
in turn, substantially smaller than the exponent within either Hartree or Hartrcc-Fock theory (C^ ~ 0.0019 a.u.) 



25| . The DMC energy depends much less sensitively than the VMC energy on the value of C 



24|. However, the 



energy difference between DMC results obtained with C = Cdmc ^-'^'^ ^ ~ ^vmc significant, while the fixed- node 
error that results from using C = is very substantial. We find that the DMC-optimized exponents of 64-electron 
crystals at different densities are given by C^^,j = 0.071 r,, '^^^ and that the DMC-optimized exponents for 16-, 64-, 
and 196-electron crystals at = 40 a.u. are very similar. We have therefore used C = C'dmc densities and 

system sizes. For antiferromagnetic crystals we have used C ~ C'dmc ~ 0.082 rj'^^^. 

Simulations were performed with up to 162, 109, 100, and 121 electrons for the paramagnetic fluid, fully spin- 
polarized (ferromagnetic) fluid, antiferromagnetic crystal, and ferromagnetic crystal, respectively 



2JI . We eliminated 



single-particle finite-size effects from the fluid energies by twist averaging within the canonical ensemble [26| . Every so 
often during VMC or DMC simulations, an offset to the grid of G vectors was chosen at random in the first Brillouin 
zone of the simulation cell, the lowest-energy plane- wave orbitals were selected, and a short period of re-equilibration 
was carried out. The finite-size errors in the energy per electron resulting from the compression of the exchange- 
correlation hole into the simulation cell and the neglect of long-ranged correlations in the kinetic energy fall off as 



N~^^^ [l9|. We therefore extrapolated our fluid energies to infinite system size by fitting our data at each density 
to En = Eao — c/N^^'^, where En is the energy per electron of the A^-electron system and c and Erx, are fitting 
parameters 24 . This differs from the form of bias that has been incorrectly assumed in previous studies of the 2D 
HEG [6|,L7|,ll4|. 

Our DMC results for the different phases are given in TableHl Our Wigner-crystal energy data were fitted to the first 
five terms in the low-density expansion of the crystal energy {E ~ ci?'^^ + C3/2fs ^^'^ + £2^^^ + c^/2i^s + c^r^'^) 27|. 
The first term is the Madelung energy of the static lattice, while the second is the quasiharmonic zero-point phonon 
energy. The corresponding coefficients can be determined analytically: ci = —1.106103 and C3/2 = 0.814 [28|. The 
remaining three coefficients were determined by fitting to our QMC data, giving C2 = 0.113743, C5/2 = —1.184994, 
and C3 = 3.097610 for the ferromagnetic crystal and C2 = 0.2662977, C5/2 = -2.63286, and C3 6.246358 for 
the antiferromagnetic crystal. We fitted our fluid energy data to the parameterization of the correlation energy 
suggested by Rapisarda and Senatore [Eq. (28) of Ref. |7[. For the paramagnetic fluid we find uq = —0.1863052 a.u., 
ai = 6.821839, 02 = 0.155226, and 03 = 3.423013, where oq, ai, 02, and 03 are the free parameters in Rapisarda and 
Senatore's fitting form. For our ferromagnetic data we obtain qq ~ —0.2909102 a.u., ai = —0.6243836, a2 = 1.656628, 



and a3 = 3.791685. Wc find a greater difference between the energies of the paramagnetic and ferromagnetic Fermi 
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DMC energy (10^^ a.u. / elec. 





Ferro. crystal Antif. crys. 


Para, fluid 


Ferro. fluid 


15 


-5.966 5(1) 








20 


-4.619 5(2) 


—4.622 9(2) 


-4.630 5(4) 


-4.6213(3) 


25 


-3.773 1(2) 


-3.775 1(3) 


-3.7774(2) 


-3.774 0(2) 


30 


—3.191 7(2) 


—3.192 2(2) 


-3.192 6(1) 


-3.1913(1) 


35 


-2.766 9(1) 


-2.767 2(1) 


-2.766 5(1) 


-2.765 7(1) 


40 


-2.443 2(1) 


-2.4431(2) 


-2.4416(1) 


-2.4416(1) 


45 


-2.1881(1) 


-2.187 5(2) 






50 


-1.9817(2) 


-1.9814(2) 







TABLE I: DMC energy as a function of for ferromagnetic and antiferromagnetic triangular Wigner crystals, and paramagnetic 
and fully ferromagnetic Fermi fluids. All results have been extrapolated to zero time step and infinite system size. 



fluids at Vs = 20 a.u. than Attaccalite et al. IJ], resulting in somewhat smaher values of the spin susceptibility. 

The ground-state energies of the different phases of the 2D HEG are plotted against rs in Fig. [TJ Unlike previous 
QMC studies, our statistieal error bars are sufhciently small that we can resolve the energy difference between the 
ferromagnetic and paramagnetic fluids. We can also identify the crystallization density with much greater precision. 



6| and = 34(4) a.u. 



31 . We find the crystallization 



Previous studies found crystallization to occur at = 37(5) a.u 
density to be = 31(1) a.u., but the transition is from the paramagnetic fiuid to an antiferromagnetic crystal, not the 
ferromagnetic fiuid to the ferromagnetic crystal as found in the previous studies. Our calculations locate the density 
at which the paramagnetic fiuid becomes unstable to the ferromagnetic fluid at about Vs = 40 a.u., but at this density 
the Wigner crystal is more stable than either fluid phase, so we do not find a region of stability for the ferromagnetic 
fiuid. We have looked for a region of stability for a partially spin-polarized fiuid by calculating the energy of a fluid 
of spin polarization ~ 2/5 at — 35 a.u. The DMC energy, extrapolated to zero time step and inflnite size, is 
—0.027666(1) a.u. per electron, which is not significantly less than the energies of the paramagnetic and ferromagnetic 
fiuids, but is significantly higher than the crystal energy (see Table|l|. It is therefore unlikely that a region of stability 
exists for a partially spin-polarized fluid. We do, however, find a transition from an antiferromagnetic Wigner crystal 
to a ferromagnetic one at = 38(5) a.u. This is considerably higher than the density at which Bernu et al. [2( 
found a transition from a spin liquid to a ferromagnetic crystal using a multispin exchange model (r^ = 175 a.u. 
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FIG. 1: (Color online) DMC energy as a function of for ferromagnetic and antiferromagnetic triangular Wigner crystals 
and paramagnetic and fully ferromagnetic Fermi fluids. In each case the Madelung energy of a triangular crystal has been 
subtracted from the energy and the result has been rescaled by rl^^ to highlight the differences between the curves. 



29[, which is somewhat lower than the 



An experimental result for the crystallization density is Vs = 35.1(9) a.u. 
QMC crystallization density. This suggests that the ideal 2D HEG is not a perfect model for electron layers in real 
semiconductor devices. 

It has been argued that the transition from a 2D Fermi fluid to a Wigner crystal cannot be first order, because, at 
the transition density, it is energetically favorable to create boundaries between macroscopically separated regions of 
fluid and crystal 8| . The energy of the Bose fluid is substantially lower than that of the Fermi fluid in the vicinity 



of the crystallization density in 2D, unlike 3D 



30[. There is therefore more scope for interesting phase behavior in 



2D. Various intermediate phases have been proposed in the literature, such as a hexatic phase with orientational 



but not translational order [9|, a supersolid phase 



10[, or a microemulsion phase [8[. Falakshahi and Waintal 
have suggested that a "hybrid" phase is stable in the vicinity of the transition density from a ferromagnetic fluid 
to a ferromagnetic Wigner crystal. (We have found that a paramagnetic fluid is stable at this density; however, we 



restrict our attention to ferromagnetic HEGs to investigate the proposed hybrid phase.) The hybrid phase has the 
same symmetry as the Wigner crystal, but has partially delocalized orbitals. 

Falakshahi and Waintal generated hybrid orbitals for the ferromagnetic TV-electron system by solving the 
Schrodinger equation for a single, positively charged particle moving in a lattice of negative charges placed at the N 
Wigner-crystal lattice sites. If the charge of the test particle is zero then the orbitals for the system are plane waves, 
i.e., fluid orbitals are obtained. If the charge is large then localized (crystal-like) orbitals are obtained. As the charge 
of the test particle is increased from zero, there must come a point at which a band gap opens up between the N and 
{N + l)th states. The set of orbitals obtained at the point at which the gap opens correspond to the hybrid phase. 
We have tried to find such a hybrid phase by optimizing orbitals of the form 

p GeP 

centered on the crystal lattice sites, where P denotes a star of symmetry-equivalent simulation-cell G vectors. The 
coefHcients {cp} and C are optimizable parameters. Equation ([T]) is a general expansion of Wannier-like orbitals. It 
can therefore describe the crystal and hybrid phases, but not the fluid, since the latter corresponds to a partially filled 
band. 

Our QMC results for ferromagnetic HEGs at = 30 a.u. obtained with different orbitals are shown in Table 
HI Within VMC it is possible to lower the energy significantly by optimizing the plane- wave coefficients in Eq. (fTj). 
apparently suggesting that Gaussian crystal orbitals are nonoptimal. This does not carry over to DMC, however. 
In fact it is more important to use Gaussian exponents optimized within DMC than it is to use either plane-wave 
expansions in the orbitals or BE functions, suggesting that the fixed-node errors in our crystal DMC energies are very 
small. We have searched for the hybrid phase by optimizing wave functions with orbitals of the form given in Eq. ([T]) 
using different starting points, but have not been able to lower the DMC energy significantly. While this does not 
prove the nonexistence of the hybrid phase, as searching for minima in a high-dimensional space is a very difficult 
problem, the fact that our extensive searches have been unable to find it makes its existence unlikely. 

In summary, we have studied the zero-temperature phase behavior of the 2D HEG using QMC. We find a transition 
from a paramagnetic fiuid to an antiferromagnetic triangular Wigner crystal at = 31(1) a.u. We find no region of 
stability for a ferromagnetic fluid, although we flnd that the paramagnetic fluid is unstable to the ferromagnetic fiuid 
at about rs = 40 a.u. and the antiferromagnetic crystal is unstable to the ferromagnetic crystal at = 38(5) a.u. 



We find no evidence for the existence of hybrid phases of the type suggested by Falakshahi and Waintal 



ll| , but we 
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Method 


Orbitals 


BF 


Energy (a.u. / elec 


VMC 


DMC-opt. 


Gauss. 


No 


-0.031849 8(1) 


VMC 


DMC-opt. Gauss.+PW 


No 


-0.031849 9(1) 


VMC 


VMC-opt. 


Gauss. 


No 


-0.031856 2(1) 


VMC 


Gauss. - 


-PW 


No 


-0.0318619(3) 


VMC 


Gauss.- 


-PW 


Yes 


-0.0318713(1) 


DMC 


Gauss.- 


-PW 


Yes 


-0.031918 0(3) 


DMC 


Gauss.- 


-PW 


No 


-0.031918 4(3) 


DMC 


VMC-opt. 


Gauss. 


No 


-0.031918 4(3) 


DMC 


DMC-opt. 


Gauss. 


No 


-0.031919 0(2) 


DMC 


DMC-opt. Gauss.+PW 


No 


-0.0319201(3) 



TABLE II: Non-twist-averaged QMC energies for a ferromagnetic, 121-electron HEG at rs = 30 a.u. obtained using different 
orbitals: crystal (Gaussian) orbitals, in which the Gaussian exponent has been optimized within VMC or DMC and "hybrid" 
(Gauss. H-plane-wave) orbitals of the form given in Eq. ([T]) with the coefficients of 20 stars of G vectors optimized within VMC. 
The orbitals were centered on the lattice sites of a triangular crystal. VMC-optimized BF functions were used in some of the 
calculations. DMC energies were obtained using a time step of 1 a.u. The non-twist-averaged DMC energy of the ferromagnetic 
Fermi fluid at this system size and density [—0.0319354(5) a.u. per electron] is relatively low due to single-particle flnite-size 
effects, so that, if anything, flnite-size effects are expected to favor the hybrid phase. 

cannot, of course, rule out the existence of other types of intermediate phase. 
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